
#ifndef _LINEAR_H
#define _LINEAR_H_

#include "Tuna.hpp"

namespace Tuna {

  template <typename T, int Dim> class Linear;

  template <typename T>
  struct Typeinfo <Linear <T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

  template <typename T>
  struct Typeinfo <Linear <T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };

template <typename Tprec, int Dim>
class Linear 
{
///
///  Traits to define the type of arrays to be used.
///
///  floatTinyArray_t : Array of floats mostly used for deltas of meshes
///
//    typedef Tprec prec_t;
    typedef typename TunaArray<Tprec, Dim>::tiny floatTinyArray_t;
    typedef typename TunaArray<Tprec, Dim>::huge ScalarField;
    
 public:
    Linear() { };
    ~Linear() { };

/*!
 * This function locate the position of the particle on the 2D mesh. 
 * The cell where the particle is located is identified with the indexes 
 * (I, J). The local position inside the cell (I,J) is defined by 
 * (alpha, beta).
 * \param pos Position of the particle.
 * \param deltas Deltas on each axis of the mesh.
 */
    inline void locate2D(floatTinyArray_t &pos, const floatTinyArray_t &deltas)
	{
	    Tprec a = pos(0) / deltas (0);
	    Tprec b = pos(1) / deltas (1);
	    
	    I = static_cast<int>(a);
	    J = static_cast<int>(b);
	    
	    alpha = a - I;
	    beta  = b - J;
	}

    inline void locate3D(floatTinyArray_t &pos, const floatTinyArray_t &deltas)
	{
	    Tprec a = pos(0) / deltas (0);
	    Tprec b = pos(1) / deltas (1);
	    Tprec c = pos(2) / deltas (2);
	    
	    I = static_cast<int>(a);
	    J = static_cast<int>(b);
	    K = static_cast<int>(c);
	    
	    alpha = a - I;
	    beta  = b - J;
	    gamma = c - K;
	}

/*!
 * This is a bilinear interpolation of the velocity to the actual position
 * of the particle. The result is stored in the array UVW.
 * \param pos Position of the particle.
 * \param deltas Deltas on each axis of the mesh.
 * \param u X-component of the velocity.
 * \param v Y-component of the velocity.
 */   
    inline floatTinyArray_t interpolate(floatTinyArray_t &pos,
					const floatTinyArray_t &deltas,
					ScalarField &u,
					ScalarField &v)
	{
	    locate2D(pos, deltas);
	    
	    UVW(0) = ( u(I  ,J  ) * (1.0 - alpha) +
		       u(I+1,J  ) * alpha ) * (1.0 - beta) +
		     ( u(I  ,J+1) * (1.0 - alpha) +
		       u(I+1,J+1) * alpha ) * beta;
	    
	    UVW(1) = ( v(I  ,J  ) * (1.0 - alpha) +
		       v(I+1,J  ) * alpha ) * (1.0 - beta) +
		     ( v(I  ,J+1) * (1.0 - alpha) +
		       v(I+1,J+1) * alpha ) * beta;
	    return UVW;
	}

    inline floatTinyArray_t interpolate(floatTinyArray_t &pos,
					const floatTinyArray_t &deltas,
					ScalarField &u,
					ScalarField &v,
					ScalarField &w)
	{
	    locate3D(pos, deltas);

	    UVW(0) = ( ( u(I  ,J  ,K  ) * (1 - alpha) +
			 u(I+1,J  ,K  ) * alpha ) * (1 - beta) +
		       ( u(I  ,J+1,K  ) * (1 - alpha) +
			 u(I+1,J+1,K  ) * alpha ) * beta ) * (1 - gamma) +
		     ( ( u(I  ,J  ,K+1) * (1 - alpha) +
			 u(I+1,J  ,K+1) * alpha ) * (1 - beta) +
		       ( u(I  ,J+1,K+1) * (1 - alpha ) +
			 u(I+1,J+1,K+1) * alpha ) * beta ) * gamma;
	    
	    UVW(1) = ( ( v(I  ,J  ,K  ) * (1 - alpha) +
			 v(I+1,J  ,K  ) * alpha ) * (1 - beta) +
		       ( v(I  ,J+1,K  ) * (1 - alpha) +
			 v(I+1,J+1,K  ) * alpha ) * beta ) * (1 - gamma) +
		     ( ( v(I  ,J  ,K+1) * (1 - alpha) +
			 v(I+1,J  ,K+1) * alpha ) * (1 - beta) +
		       ( v(I  ,J+1,K+1) * (1 - alpha ) +
			 v(I+1,J+1,K+1) * alpha ) * beta ) * gamma;
		     
	    UVW(2) = ( ( w(I  ,J  ,K  ) * (1 - alpha) +
			 w(I+1,J  ,K  ) * alpha ) * (1 - beta) +
		       ( w(I  ,J+1,K  ) * (1 - alpha) +
			 w(I+1,J+1,K  ) * alpha ) * beta ) * (1 - gamma) +
		     ( ( w(I  ,J  ,K+1) * (1 - alpha) +
			 w(I+1,J  ,K+1) * alpha ) * (1 - beta) +
		       ( w(I  ,J+1,K+1) * (1 - alpha ) +
			 w(I+1,J+1,K+1) * alpha ) * beta ) * gamma;
	    
	    return UVW;
	}

 private:
  int I, J, K;
  Tprec alpha, beta, gamma;
  floatTinyArray_t UVW;
};

} // Tuna namespace

#endif // _LINEAR_H_
